rm(list=ls())

### Loading Functions
source("functions/pica2.R")
source("functions/pica1_use.R")
source("functions/bounds_frechet_function_for_pica2.R")

# install the spsR package from naoki-egami/spsR. For reproducibility, we specify the exact version.
library(devtools)
install_github("naoki-egami/spsR@1532c8b5445ed282f203410bbad76345c69d71bf", 
               dependencies = TRUE)
library(spsR)
library(tidyverse)
library(modelsummary)
library(stargazer)

## ##########################################
## Appendix F: External Validity Analysis 
## ##########################################
# load data
load(file = "data/pica2_data.rdata")

uniq_country <- sort(unique(data$country))
pol_est_mat_1 <- pol_est_mat_2 <- matrix(NA, nrow = 5, ncol = 4)
econ_est_mat_1 <- econ_est_mat_2 <- matrix(NA, nrow = 5, ncol = 4)

for(i in 1:5){
  data_sub <- subset(data, country == uniq_country[i])
  data_sub_use <- data_sub[data_sub$D == 1, ]
  data_sub_use$A_use <- factor(data_sub_use$A, levels = c(3, 1, 2))
  
  lm_fit_pol <- lm(political_model ~ A_use, data = data_sub_use)
  pol_est_mat_1[i, 1:4] <- summary(lm_fit_pol)$coef[2, ]
  pol_est_mat_2[i, 1:4] <- summary(lm_fit_pol)$coef[3, ]
  
  lm_fit_econ <- lm(econ_model ~ A_use, data = data_sub_use)
  econ_est_mat_1[i, 1:4] <- summary(lm_fit_econ)$coef[2, ]
  econ_est_mat_2[i, 1:4] <- summary(lm_fit_econ)$coef[3, ]
  
}
rownames(pol_est_mat_1) <- rownames(pol_est_mat_2) <- rownames(econ_est_mat_1) <- rownames(econ_est_mat_2) <- uniq_country

### Figure F2 
pdf("figures/Figure_F2.pdf", height = 8, width = 10)
par(mfrow = c(2,2))
plot(seq(1:5), pol_est_mat_1[, 1], pch = 19, ylim = c(-3, 3), main = "Effects of US Propaganda \n on Political Models",
     xaxt = "n", xlab = "Sites", ylab = "Estimates")
segments(seq(1:5), pol_est_mat_1[, 1] - 1.96*pol_est_mat_1[, 2], 
         seq(1:5), pol_est_mat_1[, 1] + 1.96*pol_est_mat_1[, 2])
abline(h = 0, lty = 2, lwd = 2)
Axis(side = 1, at = seq(1:5), labels = uniq_country)

plot(seq(1:5), econ_est_mat_1[, 1], pch = 19, ylim = c(-3, 3), main = "Effects of US Propaganda \n on Economic Models",
     xaxt = "n", xlab = "Sites", ylab = "Estimates")
segments(seq(1:5), econ_est_mat_1[, 1] - 1.96*econ_est_mat_1[, 2], 
         seq(1:5), econ_est_mat_1[, 1] + 1.96*econ_est_mat_1[, 2])
abline(h = 0, lty = 2, lwd = 2)
Axis(side = 1, at = seq(1:5), labels = uniq_country)

plot(seq(1:5), pol_est_mat_2[, 1], pch = 19, ylim = c(-3, 3), main = "Effects of Chinese Propaganda \n on Political Models",
     xaxt = "n", xlab = "Sites", ylab = "Estimates")
segments(seq(1:5), pol_est_mat_2[, 1] - 1.96*pol_est_mat_2[, 2], 
         seq(1:5), pol_est_mat_2[, 1] + 1.96*pol_est_mat_2[, 2])
abline(h = 0, lty = 2, lwd = 2)
Axis(side = 1, at = seq(1:5), labels = uniq_country)

plot(seq(1:5), econ_est_mat_2[, 1], pch = 19, ylim = c(-3, 3), main = "Effects of Chinese Propaganda \n on Economic Models",
     xaxt = "n", xlab = "Sites", ylab = "Estimates")
segments(seq(1:5), econ_est_mat_2[, 1] - 1.96*econ_est_mat_2[, 2], 
         seq(1:5), econ_est_mat_2[, 1] + 1.96*econ_est_mat_2[, 2])
abline(h = 0, lty = 2, lwd = 2)
Axis(side = 1, at = seq(1:5), labels = uniq_country)
dev.off()

### SPS Estimator 
# First, we import the SPS site selection results. 
load(file = "sps_out_pica2_final.rdata")

sps_est_pol_US <- sps_estimator(out = out, estimates_selected = pol_est_mat_1[, c(1,2)])
sps_est_econ_US <- sps_estimator(out = out, estimates_selected = econ_est_mat_1[, c(1,2)])
sps_est_pol_China <- sps_estimator(out = out, estimates_selected = pol_est_mat_2[, c(1,2)])
sps_est_econ_China <- sps_estimator(out = out, estimates_selected = econ_est_mat_2[, c(1,2)])

sps_est_mat <- rbind(sps_est_pol_US$average_site_ATE, sps_est_econ_US$average_site_ATE, 
                     sps_est_pol_China$average_site_ATE, sps_est_econ_China$average_site_ATE)

id_use <- c(1.75, 2.25, 2.75, 3.25)

pdf("figures/Figure_F3.pdf", height = 5, width = 6)
plot(id_use, sps_est_mat[, 1], pch = 19, ylim = c(-3, 4), xlim = c(1.5, 3.5), 
     main = "Average-Site ATEs",
     xaxt = "n", xlab = "Outcomes", ylab = "Estimates")
segments(id_use, sps_est_mat[, 1] - 1.96*sps_est_mat[, 2], 
         id_use, sps_est_mat[, 1] + 1.96*sps_est_mat[, 2])
abline(h = 0, lty = 2, lwd = 2)
Axis(side = 1, at = id_use, labels = c("Political\n Models", "Economic\n Models",
                                       "Political\n Models", "Economic\n Models"), line = 0.5, tick = 0, lwd = 0)
text("Effects of \n US Propaganda", x = 2, y = 3.5, font = 2)
text("Effects of \n Chinese Propaganda", x = 3, y = 3.5, font = 2)
dev.off()

### Site-level Cross-Validation
# This helps assess the influence of unobserved confounders. 
cv_out_pol_1 <- sps_cv(out = out, estimates_selected = pol_est_mat_1[,c(1,2)])

cv_out_econ_1 <- sps_cv(out = out, estimates_selected = econ_est_mat_1[,c(1,2)])

cv_out_pol_2 <- sps_cv(out = out, estimates_selected = pol_est_mat_2[,c(1,2)])

cv_out_econ_2 <- sps_cv(out = out, estimates_selected = econ_est_mat_2[,c(1,2)])

round(c(cv_out_pol_1$p_value, cv_out_econ_1$p_value), 2)
round(c(cv_out_pol_2$p_value, cv_out_econ_2$p_value), 2)


## ##########################
## Appendix G
## ##########################

##
## Attrition
##

#Incompletes
load(file = "data/pica2_incompletes.rdata")

#Attrition by group
lm1 <- lm_robust(attrited_all~as.factor(Group), data=d)

lm2 <- lm_robust(attrited_all~Treatment, data=d, subset = 
                   Treatment=="Forced Placebo"|Treatment=="China Forced Treatment")

lm3 <- lm_robust(attrited_all~Treatment, data=d, subset = 
                   Treatment=="Forced Placebo"|Treatment=="US Forced Treatment")

lm4 <- lm_robust(attrited_all~Treatment, data=d, subset = 
                   Treatment=="China Choice Placebo"|Treatment=="China Choice Treatment")

lm5 <- lm_robust(attrited_all~Treatment, data=d, subset = 
                   Treatment=="US Choice Placebo"|Treatment=="US Choice Treatment")

lm6 <- lm_robust(attrited_all~Treatment, data=d, subset = 
                   Treatment=="Free Nature"|Treatment=="US Choice Treatment"|Treatment=="China Choice Treatment")


# Summarize attrition
modelsummary(list("(1)"=lm1, "(2)"=lm2, "(3)"=lm3, "(4)"=lm4, "(5)"=lm5),
             coef_map=c("as.factor(Group)Forced"="Forced v. Choice",
                        "TreatmentChina Forced Treatment"="China Forced v. Forced Placebo",
                        "TreatmentUS Forced Treatment"="US Forced v. Forced Placebo",
                        "TreatmentChina Choice Treatment"="China Choice v. China Choice Placebo",
                        "TreatmentUS Choice Treatment"="US Choice v. US Choice Placebo",
                        "(Intercept)"="Intercept"),
             output = "figures/Table_G1.tex")

rm(d)

load("data/pica2_data_include_ac_failure.rdata")

###
### Attention Check
###
m1 <- lm(both_passed ~ female + education + age + ideology_right + national_pride, data=data)
m2 <- lm(attentioncheck1_pass ~ female + education + age + ideology_right + national_pride, data=data)
m3 <- lm(attentioncheck2_pass~ female + education + age + ideology_right + national_pride, data=data)

stargazer(m1, m2, m3, type = "latex", out = "figures/Table_G2.tex")

###
### Analysis including people who failed attention checks
###

### Political Model as the Outcome Variable 
out_political_F <- pica2(Y = "political_model", D = "D", A = "A", C = "actual_choice", 
                       covariates = c("female", "education", "national_pride", "ideology_right", "country"),
                       placebo = "placebo", placebo_level = 3, data = data)
rownames(out_political_F$ATE) <- c("US vs Nature", "China vs Nature")  

rownames(out_political_F$ACTE_view) <- c("Effect of US video on US-viewer", "Effect of China video on China-viewer")
rownames(out_political_F$ACTE_non_view) <- c("Effect of US video on non-US-viewer", "Effect of China video on non-China-viewer")

### Economic Model as the Outcome Variable 
out_econ_F <- pica2(Y = "econ_model", D = "D", A = "A", C = "actual_choice", 
                  covariates = c("female", "education", "national_pride", "ideology_right", "country"),
                  placebo = "placebo", placebo_level = 3, data = data)
rownames(out_econ_F$ATE) <- c("US vs Nature", "China vs Nature")  

rownames(out_econ_F$ACTE_view) <- c("Effect of US video on US-viewer", "Effect of China video on China-viewer")
rownames(out_econ_F$ACTE_non_view) <- c("Effect of US video on non-US-viewer", "Effect of China video on non-China-viewer")

### World Leader as the Outcome Variable 
out_world_F <- pica2(Y = "world_leader", D = "D", A = "A", C = "actual_choice", 
                   covariates = c("female", "education", "national_pride", "ideology_right", "country"),
                   placebo = "placebo", placebo_level = 3, data = data)
rownames(out_world_F$ATE) <- c("US vs Nature", "China vs Nature")  

rownames(out_world_F$ACTE_view) <- c("Effect of US video on US-viewer", "Effect of China video on China-viewer")
rownames(out_world_F$ACTE_non_view) <- c("Effect of US video on non-US-viewer", "Effect of China video on non-China-viewer")


pt <- as.data.frame(rbind(out_political_F$ATE, out_political_F$ACTE_view, out_political_F$ACTE_non_view))
et <- as.data.frame(rbind(out_econ_F$ATE, out_econ_F$ACTE_view, out_econ_F$ACTE_non_view))
wt <- as.data.frame(rbind(out_world_F$ATE, out_world_F$ACTE_view, out_world_F$ACTE_non_view))

pt <- pt[, c("Estimate", "Pr(>|t|)")]
et <- et[, c("Estimate", "Pr(>|t|)")]
wt <- wt[, c("Estimate", "Pr(>|t|)")]

t <- cbind(pt, et, wt)
library(xtable)
t <- xtable(t)

# Table G3
print(t, file = "figures/Table_G3.tex", type = "latex")


## manipulation check 


load("data/pica2_data_include_ac_failure.rdata")



# Ensure 'mc' is a character column
data <- data %>% mutate(mc = as.character(mc))

# Create dummy variables for each response category
data <- data %>%
  mutate(
    mc_problem  = as.integer(grepl("problem", mc, ignore.case = TRUE)),
    mc_nature   = as.integer(grepl("nature", mc, ignore.case = TRUE)),
    mc_american = as.integer(grepl("American society", mc, ignore.case = TRUE)),
    mc_chinese  = as.integer(grepl("Chinese society", mc, ignore.case = TRUE)),
    mc_japanese = as.integer(grepl("Japanese society", mc, ignore.case = TRUE))
  )

# Check if the dummies were created correctly
table(data$mc_problem)
table(data$mc_nature)
table(data$mc_american)
table(data$mc_chinese)
table(data$mc_japanese)


# Define experimental conditions to analyze
treatment_conditions <- c("china treatment_forced", "china treatment_free", "china nature_free",
                          "us treatment_forced", "us treatment_free", "us nature_free",
                          "nature_forced", "nature_free")

# Function to calculate manipulation check proportions for each treatment condition
calculate_mc_proportions <- function(data, treatment_condition, mc_vars) {
  category_data <- data %>% filter(treatment == treatment_condition)
  
  # Compute proportions
  result <- category_data %>%
    summarize(across(all_of(mc_vars), ~ mean(.x, na.rm = TRUE))) %>%
    pivot_longer(cols = everything(), names_to = "dummy_variable", values_to = "prop_result") %>%
    mutate(treatment = treatment_condition)
  
  return(result)
}

# Apply function to all treatment conditions
mc_results <- bind_rows(lapply(treatment_conditions, calculate_mc_proportions, data = data, 
                               mc_vars = c("mc_problem", "mc_nature", "mc_american", "mc_chinese", "mc_japanese")))

# Reshape data to create the LaTeX table format
mc_table <- mc_results %>%
  pivot_wider(names_from = treatment, values_from = prop_result) %>%
  rename(Type = dummy_variable)

# Convert to LaTeX format
tbl <- xtable(mc_table)

# Export the table to a LaTeX file
output_file <- "mc.tex"
print(tbl, file = output_file, include.rownames = FALSE, floating = FALSE)

print(tbl, file = "figures/Table_G4.tex", type = "latex")


## ##########################################
## Appendix H: Full Results
## ##########################################
# reload clean data
load(file = "data/pica2_data.rdata")

### ############
### Section H.1
### ############
out_political <- pica2(Y = "political_model", D = "D", A = "A", C = "actual_choice", 
                       covariates = c("female", "education", "national_pride", "ideology_right", "country"),
                       placebo = "placebo", placebo_level = 3, data = data)
out_econ <- pica2(Y = "econ_model", D = "D", A = "A", C = "actual_choice", 
                  covariates = c("female", "education", "national_pride", "ideology_right", "country"),
                  placebo = "placebo", placebo_level = 3, data = data)
out_world <- pica2(Y = "world_leader", D = "D", A = "A", C = "actual_choice", 
                   covariates = c("female", "education", "national_pride", "ideology_right", "country"),
                   placebo = "placebo", placebo_level = 3, data = data)


## combining all results to export full tables for appendix 
pt <- as.data.frame(rbind(out_political$ATE, out_political$ACTE_view, out_political$ACTE_non_view))
et <- as.data.frame(rbind(out_econ$ATE, out_econ$ACTE_view, out_econ$ACTE_non_view))
wt <- as.data.frame(rbind(out_world$ATE, out_world$ACTE_view, out_world$ACTE_non_view))

pt <- pt[, c("Estimate", "Pr(>|t|)")]
et <- et[, c("Estimate", "Pr(>|t|)")]
wt <- wt[, c("Estimate", "Pr(>|t|)")]

t <- cbind(pt, et, wt)
rownames(t) <- c("US vs Nature", "China vs Nature", 
                 "Effect of US video on US-viewer",
                 "Effect of China video on China-viewer", 
                 "Effect of US video on non-US-viewer",
                 "Effect of China video on non-China-viewer"
)
library(xtable)
t <- xtable(t)

### Table H1
print(t, file = "figures/Table_H1.tex", type = "latex")


### Table H2: Report coefficients
library(stargazer)
H2 <- stargazer(list(out_political$lm_ate_out, out_econ$lm_ate_out, out_world$lm_ate_out))
writeLines(capture.output(stargazer(H2)), "figures/Table_H2.tex")

### Figure H1: Distribution of Outcomes
tab2 <- table(data$political_model[data$actual_choice == 2 & data$D == 0 & data$placebo == 0])
tab2 <- c(0, tab2)
names(tab2)[1] <- "1"

pdf("figures/Figure_H1.pdf", height = 5, width = 9)
par(mfrow = c(1,2), oma = c(0, 0, 3, 0))
barplot(table(data$political_model[data$actual_choice == 2 & data$D == 0 & data$placebo == 1]), 
        main = "Placebo-condition", ylim = c(0, 20))
barplot(tab2, main = "Treatment-condition", ylim = c(0, 20))
mtext(side = 3, at = 0.5, outer = TRUE, text = "Political Models: China-Choosers", font = 2, cex = 1.35)
dev.off()

### ############
### Section H.2
### ############

### Table H3
library(xtable)
H3 <- xtable(round(prop.table(table(data$country,  data$stated_preference), 1), 2))
print(H3, file = "figures/Table_H3.tex", type = "latex")

### Table H4
H4 <- xtable(round(prop.table(table(data$country[data$group == "free"], 
                              data$freechoice_mediapref[data$group == "free"]), 1), 2))
print(H4, file = "figures/Table_H4.tex", type = "latex")

### Table H5
H5 <- xtable(addmargins(table(data$freechoice_mediapref[data$group == "free"],
                        data$stated_preference[data$group == "free"])), digits = 0)
print(H5, file = "figures/Table_H5.tex", type = "latex")


### ############
### Section H.3
### ############
### Subgroup Analyses
ACTE_political <- rbind(out_political$ACTE_view, out_political$ACTE_non_view)
out_political_boot <- pica2_boot(Y = "political_model", D = "D", A = "A", C = "actual_choice", 
                                 covariates = c("female", "education", "national_pride", "ideology_right", "country"),
                                 placebo = "placebo", placebo_level = 3, data = data, 
                                 boot_num =  1000, seed = 1234)

diff_US_boot <- out_political_boot$ACTE_view_boot[1, ] - out_political_boot$ACTE_non_view_boot[1, ]
diff_China_boot <- out_political_boot$ACTE_view_boot[2, ] - out_political_boot$ACTE_non_view_boot[2, ]

US_subgroup <- c(mean(diff_US_boot), sd(diff_US_boot))
China_subgroup <- c(mean(diff_China_boot), sd(diff_China_boot))

US_subgroup <- c(ACTE_political[1,1] - ACTE_political[3,1], sd(diff_US_boot))
China_subgroup <- c(ACTE_political[2,1] - ACTE_political[4,1], sd(diff_China_boot))

### Subgroup Analyses with Other Variables
# Looking at the forced group 
data_forced <- data[data$D == 1, ]
data_forced$A <- factor(data_forced$A, levels = c(3, 1, 2))

subgroup_var <- c("female", "education", "national_pride", "ideology_right", "age")
out_ATE_pol_1 <- lm_robust(political_model ~ A*female, data = data_forced)
summary(out_ATE_pol_1)

data_forced$edu_bin <- case_when(data_forced$education >= median(data_forced$education, na.rm = TRUE) ~ 1,
                                 data_forced$education < median(data_forced$education, na.rm = TRUE) ~ 0)
out_ATE_pol_2 <- lm_robust(political_model ~ A*edu_bin, data = data_forced)
summary(out_ATE_pol_2)

data_forced$np_bin <- case_when(data_forced$national_pride >= median(data_forced$national_pride, na.rm = TRUE) ~ 1,
                                data_forced$national_pride < median(data_forced$national_pride, na.rm = TRUE) ~ 0)
out_ATE_pol_3 <- lm_robust(political_model ~ A*np_bin, data = data_forced)
summary(out_ATE_pol_3)

data_forced$ideo_bin <- case_when(data_forced$ideology_right >= median(data_forced$ideology_right, na.rm = TRUE) ~ 1,
                                  data_forced$ideology_right < median(data_forced$ideology_right, na.rm = TRUE) ~ 0)
out_ATE_pol_4 <- lm_robust(political_model ~ A*ideo_bin, data = data_forced)
summary(out_ATE_pol_4)

data_forced$age_bin <- case_when(data_forced$age >= median(data_forced$age, na.rm = TRUE) ~ 1,
                                 data_forced$age < median(data_forced$age, na.rm = TRUE) ~ 0)
out_ATE_pol_5 <- lm_robust(political_model ~ A*age_bin, data = data_forced)
summary(out_ATE_pol_5)

US_effect_diff <- rbind(summary(out_ATE_pol_1)$coef[5, 1:2], 
                        summary(out_ATE_pol_2)$coef[5, 1:2], 
                        summary(out_ATE_pol_3)$coef[5, 1:2],
                        summary(out_ATE_pol_4)$coef[5, 1:2], 
                        summary(out_ATE_pol_5)$coef[5, 1:2])
China_effect_diff <- rbind(summary(out_ATE_pol_1)$coef[6, 1:2], 
                           summary(out_ATE_pol_2)$coef[6, 1:2], 
                           summary(out_ATE_pol_3)$coef[6, 1:2],
                           summary(out_ATE_pol_4)$coef[6, 1:2], 
                           summary(out_ATE_pol_5)$coef[6, 1:2])
rownames(US_effect_diff) <- subgroup_var
rownames(China_effect_diff) <- subgroup_var

US_effect_diff_p <- rbind("Media Preference" = US_subgroup, 
                          US_effect_diff)
China_effect_diff_p <- rbind("Media Preference" = China_subgroup,
                             China_effect_diff)

## Table H6
H6 <- xtable(rbind(US_effect_diff_p, China_effect_diff_p))
print(H6, file = "figures/Table_H6.tex", type = "latex")
